﻿using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace lesson1
{
    class CubicSpline
    {
        public List<Point> points;
        Vector<double> _m;
        public CubicSpline(List<Point> ps)
        {
            points = ps;
            computeM();
        }

        void computeM()
        {
            int N = points.Count;
            var A = Matrix<double>.Build.Dense(N, N,0);
            var b = Vector<double>.Build.Dense(N);
           
            for (int i = 1; i < N-1; i++)
            {
                
                A[i, i - 1] = lamd(i);
                A[i, i] = 2;
                A[i, i + 1] = mu(i);
                b[i] = R(i);
            }
            A[0, 0] = 1;
            A[0, 1] = 1;
            b[0] = 2*(y(2)-y(1))/h(1);
            A[N-1, N-2] =1;
            A[N - 1, N - 1] = 1;
            b[N-1] =2*(y(N-1)-y(N-2))/h(N-2);

            _m = A.Solve(b);
        
        }
        public void draw(IntPtr hDC)
        {
            MyGDIClass.BresenhamCircle(hDC, (int)x(points.Count - 1), (int)y(points.Count - 1), 5);
            for (int i = 0; i < points.Count-1; i++)
            {
                MyGDIClass.BresenhamCircle(hDC,(int)x(i),(int)y(i),5);
                int detX =Math.Abs( (int)(x(i+1)-x(i)));

                for (int j = 0; j < detX; j++)
                {
                    double xx = x(i) + j;
                    double x_xi = xx - x(i);
                    double y = a(i) + b(i) * x_xi + c(i) * x_xi * x_xi + d(i) * x_xi * x_xi * x_xi;
                    MyGDIClass.SetPixel(hDC, (int)xx, (int)y, 0);
                }
            }

        }

        private double d(int i)
        {
            return (m(i) + m(i + 1)) / (h(i) * h(i)) - 2 * (y(i + 1) - y(i)) / (h(i) * h(i) * h(i));
        }

        private double c(int i)
        {
            return 3 * (y(i + 1) - y(i)) / (h(i) * h(i)) - (2 * m(i) + m(i + 1)) / h(i);
        }

        private double b(int i)
        {
            return m(i);
        }

        private double m(int i)
        {
            return _m[i];
        }

        private double a(int i)
        {
            return y(i);
        }

        private double lamd(int i)
        {
            return h(i) / (h(i) + h(i - 1));
        }

        private double h(int i)
        {
            return x(i + 1) - x(i);
        }

        private double x(int i)
        {
            return points[i].X;
        }
        private double y(int i)
        {
            return points[i].Y;
        }
        private double R(int i)
        {
            return 3*(mu(i)*(y(i+1)-y(i))/h(i)+lamd(i)*(y(i)-y(i-1))/h(i-1));
        }

        private double mu(int i)
        {
            return 1 - lamd(i);
        }
    }
}
